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Time-resolved hadronic particle acceleration in the 
recurrent nova RS Ophiuchi 


H.E.S.S. Collaboration? 


Recurrent novae are repeating thermonuclear explosions in the outer layers 
of white dwarfs, due to the accretion of fresh material from a binary compan- 
ion. The shock generated by ejected material slamming into the companion 
star's wind, accelerates particles to very-high-energies. We report very-high- 
energy (VHE, = 100 GeV) gamma rays from the recurrent nova RS Ophiuchi 
up to a month after its 2021 outburst, using the High Energy Stereoscopic 
System. The VHE emission has a similar temporal profile to lower-energy 
GeV emission, indicating a common origin, with a two-day delay in peak flux. 
These observations constrain models of time-dependent particle energization, 
favouring a hadronic emission scenario over the leptonic alternative. This con- 
firms that shocks in dense winds provide favourable environments for efficient 


cosmic-ray acceleration to very-high-energies. 
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RS Ophiuchi (RS Oph) is a recurrent nova system comprising a white dwarf and a com- 
panion red giant star. Novae are a source of high-energy particles (/, 2), with non-thermal 
gamma-ray emission in the range ~100 MeV to ~10 GeV (3). The RS Oph system is located 
approximately 1.4kpc from Earth (4); analysis of Gaia data suggests larger distances of 2.3 
kpc (5) or 2.7 kpc (6), although the reliability of these estimates is questionable due to the or- 
bital motion in RS Oph. The binary components have a separation of 1.48 astronomical units 
(au), close enough for the white dwarf to continually accrete material from its companion (7). 
At irregular intervals, enough material accumulates on the surface of the white dwarf to trigger 
a thermonuclear explosion, driving a quasi-spherical shock into the red giant's wind (see 8, Fig- 
ure 2). Eight outbursts were observed between 1898 and 2006, recurring in intervals of 9 to 26 
years (9). 

On Si August 2021, an outburst of RS Oph was identified from optical observations (10), 
with a peak naked-eye visual magnitude of 4.5, over a thousand times brighter than the quiescent 
visual magnitude of 12.5. We observed RS Oph with the High Energy Stereoscopic System 
(H.E.S.S.), an array of five atmospheric Cherenkov telescopes. Observations commenced on 9°? 
August 2021 and continued for five nights until 13° August 2021. Optical background emission 
from the Moon prevented good quality observations for the following ten nights. During each 
of those five nights, H.E.S.S. detected point-like gamma-ray emission from the direction of 
RS Oph, with a significance of > 6 sigma on each night (see 8, Table S1). The data for those 
five nights combined are shown in Figure |I] Observations recommenced on 25'^ August 2021, 
17 days after the initial outburst. We find evidence for a much weaker signal ( 3c above the 
background) was seen in ~15 hours (after quality cuts) of data accumulated over the following 
14 days. 

We performed a spectral analysis of the H.E.S.S. data for the first five observation nights 


separately. We also separated the data from the array of four 106 m? mirror area telescopes 
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Figure 1: RS Oph significance maps. Significance maps derived from the H.E.S.S. > 100GeV gamma-ray 
observations for the early (A) and late (B) phases of the RS Oph 2021 outburst. To = Modified Julian Day (MJD) 
59435.25, is the time of peak optical emission. The dashed white circles indicate the point-spread-function (PSF). 
(designated CT1-4) and the fifth 612 m? mirror area low-threshold telescope (CT5). We find 
that the VHE flux is variable, with a spectral index > 3 throughout (see 8, Table S2). 
Figure[2|shows the time evolution of the gamma-ray flux curve for photon energies between 
250 GeV and 2.5 TeV. The VHE gamma-ray flux rises smoothly from To, the time of peak opti- 
cal emission in the V band (/0), until a VHE peak on the third night of observations, after which 
the VHE gamma-ray energy flux decays by an order of magnitude over a two-week period. We 
obtained 60 MeV — 500 GeV data taken by the Fermi-LAT (Large Area Telescope) instrument 
for the same time period as the H.E.S.S. observations which are also shown in Figure] The 
flux varies in the range 1x 1078 — 2x10~'° erg cm ?s !, with a peak flux in the Fermi-LAT 
data on Tọ + 1 day. The VHE gamma-ray emission peak is delayed by a further two days. 
After the peak flux, we fitted the decay in time ¢ of the energy flux with a power-law with 
exponent o, t * and found best-fitting values of a ~ 1.3 — 1.4 in both data sets: aypss = 


1.43 + 0.18 for HESS and azar = 1.31 + 0.07 for Fermi-LAT, for the choice of Tọ —1 day. 
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Figure 2: Gamma-ray light curves of RS Oph. Light curves of gamma-ray emission from RS Oph including data 
from Fermi-LAT and H.E.S.S. observations. The H.E.S.S. data (red squares) cover a period of five nights, after 
which observations ceased for ten days due to bright moonlight, marked by the shaded grey band, then recom- 
menced for a period of 14 days. The H.E.S.S. flux is integrated from 250 GeV to 2.5 TeV, whilst the Fermi-LAT 
flux is integrated from 60 MeV to 500 GeV. Fermi-LAT data are shown in 6-hour bins (blue circles) corresponding 
to the time windows of the H.E.S.S. observations, and data outside of these times shown with semi-transparent 
markers. Error bars are 1 o statistical uncertainties. A power-law slope model was fitted to the temporal decay 
after the time of peak flux for both instruments (red and blue dashed lines, with uncertainties indicated by the 
shaded regions). The vertical dotted black line indicates the peak of the outburst in the optical waveband, To. 
The Fermi-LAT flux and temporal decay are consistent with that obtained from bins of 24-hour 
duration, the higher statistics of the larger bins enabling the detailed spectral analysis shown in 
Figure 

The combined H.E.S.S. and Fermi-LAT data (8) allow us to measure wide-band gamma-ray 
spectra over more than four orders of magnitude in energy and follow their temporal evolution 
(Figure 3). The RS Oph spectra are consistent with a log-parabola model. Comparison of 
spectra taken on different nights show a general trend for the flux normalisation to decrease and 
the parabola to widen over time (see 8, Table S2). 


The similarity between the spectra of the Fermi-LAT and H.E.S.S. data, and their similar 


decay profiles after their respective peaks, indicate a common origin for the gamma rays from 
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Figure 3: RS Oph gamma-ray spectra. The H.E.S.S. and Fermi-LAT spectra for 9 (green) and 13 (orange) 
August fitted with a log-parabola model. The analysis is applied separately for the H.E.S.S. CT1-4 (squares) and 
CTS (triangles) - see text. The Fermi-LAT data (open circles) are integrated over 24 h centred at the H.E.S.S. 
observation times. There is clear spectral evolution from the 9*^ to the 13*^ August, with a noticeable reduction 
in the Fermi-LAT flux as well as an increase in the maximum energy of the TeV spectrum. Error bars are 1 sigma 
statistical uncertainty, and upper limits are the 95% confidence level. 
one day to one month after the explosion. We assume that the particles that generate the gamma 
rays are accelerated at the external shock as it propagates into the wind of the red giant (8, 
Figure S2). Optical spectroscopic measurements of the 2021 nova indicate shock velocities in 
the range ug, = 4000 — 5000 kms"! (71), compatible with measurements from the previous 
2006 outburst of RS Oph (72, 13). High resolution images of the 2006 event (74) indicated the 
polar regions of the shock expanded at ~ 5000 km s^! over the first 5 months. We therefore 
assume that during the first week following the 2021 outburst the shock velocity did not fall 
below several thousand kilometers per second. 

The images of the 2006 nova showed a quasi-spherical outflow, pinched at an equatorial 


ring (/4, 15). This is consistent with a shock expanding into the wind of the red giant orthogonal 


to the orbital plane of the binary, but inhibited close to the plane by the denser gas (7, 16). We 


expect particles to undergo diffusive shock acceleration at the external fast moving shocks, 
above and below the orbital plane of the binary. We consider two scenarios to explain the 
observed spectral and temporal properties of the gamma-ray emission from RS Oph: gamma- 
ray production from accelerated protons colliding with dense gas in the downstream volume 
(hadronic 7° decay model), or gamma-ray emissions from energetic electrons scattering low 
energy photons from the nova (inverse Compton model). For both models, the observations 
place strong constraints on the physical conditions, particularly the acceleration efficiencies 
required to match the measured fluxes and maximum photon energies (8). 

VHE gamma-ray emission requires acceleration of particles to >TeV energies. The maxi- 
mum energy a particle attains at a shock is determined by either the time taken before radiative 
cooling dominates over acceleration, or by when particles become too energetic and escape up- 
stream of the shock (17). This confinement limit applies when the accelerating particles are 
unable to excite magnetic field fluctuations to a sufficient level ahead of the shock. As parti- 
cles spend longer diffusing upstream than downstream, the details of the downstream magnetic 
fields can be neglected (78). For upstream magnetic-field amplification to be effective, a suf- 
ficient flux of particles, typically protons, must escape upstream of the shock. This requires 
both efficient transfer of the shock kinetic energy to relativistic protons, and that a fraction of 
these protons penetrate far upstream. Escaping particles have energies concentrated close to the 
maximum particle energy; less energetic particles are confined to the shock. The escaping flux 


per unit area at a given shock radius can be parameterized as 
Jesc SE? [2d a) Emax b (1) 


where e is the elementary charge, F; = i Pup, is the energy flux density for a high Mach- 
number non-relativistic shock, pup is the immediate upstream gas density at that radius, Emax is 


the maximum particle energy which dominates the escaping flux, and the efficiency parameter 


Eese depends on the assumed particle spectrum (8). For a wind-like density profile, and neglect- 
ing radiative losses, the confinement limit on the maximum energy for a particle of charge |Z |e 


(with atomic number Z) is 


. 1/2 2 
Eosi M /Vwina Ush 
Emax = 1.5 Z T V 5 2 
Zi (i 101 kg m^ 5000kms 1) ^ = 


where M and Uwind are the mass-loss rate and the wind velocity of the red giant. Eese is predicted 


to be about 1% for high Mach-number shocks (17). For RS Oph, M/vwing = 6 x 10!! kg m^! 
(15) which, together with the inferred shock velocities, indicates a maximum energy Emax & 
10 TeV. This is compatible with the measured maximum photon energies Ey max ^ 1TeV 
(Figure B). 

In the hadronic scenario, the gamma-ray lightcurves are consistent with an expanding shock 
in a decreasing density profile. With the adopted distance of 1.4 kpc (4), the measured gamma- 
ray fluxes require that >10% of the post-shocked medium’s internal energy goes to accelerating 
protons or other nuclei. The delay between the peaks in the Fermi-LAT and H.E.S.S. lightcurves 
would then reflect the finite acceleration time of the >1 TeV protons, or more specifically, the 
time taken to populate the high-energy tail of the distribution (79). A simple calculation of the 
acceleration time, based on a comparison of the confinement and Hillas limits, implies that this 
should happen on the order of days (8). This is consistent with the spectral evolution seen in 
Figure B] a reduction in the Fermi-LAT flux, accompanied by a hardening in H.E.S.S. flux and 
increased E. max over the first few days post-outburst. Attenuation of gamma-rays due to the 
novas optical and infra-red photon fields is minor below 1 TeV a few hours after the explosion, 
and therefore attenuation alone cannot account for the observed hardening (see 8, Figure S10). 

For the alternative, leptonic scenario in which TeV gamma-rays are produced by VHE elec- 
trons, the acceleration needs to overcome the strong radiative losses due to inverse Compton 


cooling in the strong photon fields of the nova, as well as synchrotron cooling in the magnetic 


field in the shock region. To achieve this, electrons must accelerate at close to the Bohm rate, 1.e. 
the scattering rate equal to the rate of gyration in the magnetic field (78). Such efficient scatter- 
ing requires strong self-generated magnetic fluctuations upstream of the shock, which implies 
the presence of an energetic relativistic hadronic component. In this scenario, the differences 
between the spectral slopes in the Fermi-LAT and H.E.S.S. energy ranges are a consequence of 
the energy-dependent cooling rates in time-dependent photon fields. Electrons that radiate in 
the VHE band cool on a timescale less than the age of the nova remnant at the times the obser- 
vations were taken, while lower-energy un-cooled electrons accumulate downstream over time. 
The Fermi-LAT light curve in this scenario then reflects the evolution of the energy density of 
soft-photon targets, while the H.E.S.S. lightcurve traces the full radiative output of high-energy 
electrons up to the VHE peak. After the peak, due to the rapidly decreasing photon energy den- 
sity, the cooling time increases faster than the remnant’s age, and the VHE emitting electrons 
are also slow cooling. 

From the time-dependent numerical single-zone model, parameters can be found to approx- 
imately describe the light curves and spectra in both leptonic and hadronic scenarios, providing 
quantitative estimates for the acceleration efficiencies at early times t (t < Tọ +5 days). We are 
unable to account for the temporal decay at later times since a more complex treatment of the 
internal structure of the nova remnant, its non-spherical geometry, and escape of particles from 
the emission zone is required. Both the leptonic and hadronic models at early times are consis- 
tent with continuous injection of particles following a power law spectrum in energy E~?? with 
a high-energy cut-off. To approximately match the observed flux, the hadronic model requires 
>10% of the shocked gas’ internal energy to be transferred to non-thermal protons, while the 
leptonic model requires >1% efficiency for non-thermal electrons. 

Such a high fraction of the total energy in non-thermal electrons is inconsistent with the- 


ories of injection at high Mach-number shocks, for which ion injection efficiency is expected 


to be much higher than that of electrons (20). Numerical simulations of high Mach-number 
shocks find a ratio of «10 ? for electron to ion energy densities (27), which is consistent with 
multi-wavelength models of supernova remnants (22). A >1% efficiency of conversion to non- 
thermal electrons cannot be realised in a purely leptonic model. For this reason, we prefer the 
hadronic scenario discussed above, for which both the implied high proton acceleration efficien- 
cies and inferred maximum energy are in line with theoretical predictions (/7). Our findings 


support previous hadronic models of gamma-ray novae (23-25). 


The VHE detection of RS Oph demonstrates that particle acceleration to ~TeV energies 
can occur within the dense wind environments of recurrent novae. The total kinetic energy from 
each nova of RS Oph is estimated to be ~ 104? erg (1077 solar masses of ejecta at 4000 km s^! 
(13)), with a large fraction of this being converted to relativistic protons and heavier nuclei, 
which are the main constituents of Galactic cosmic rays. Each nova event generates enough 
cosmic rays to fill a cubic parsec (pc) volume with an energy density of ~0.1 eV cm *, similar 
to the local Galactic cosmic-ray energy density of ~ (0.8 — 1.0) eV cm ? (26) sustained by 
supernovae. In the case of RS Oph, the cosmic-ray energy input recurs approximately every 
At — 15 — 20 years, leading to an almost continuous injection of non-thermal particles. For a 
diffusion coefficient D in the neighbourhood of RS Oph, the cosmic-ray output from each nova 
is spread over a diffusion length aig = VAD At. Using a Galactic average D ~ (3 — 5) x 1078 
cm? s^! (27), we find laig > 1pc, and the contribution from each nova is sub-dominant to the 
average Galactic cosmic-ray population. If the diffusion coefficient in the neighbourhood of the 
nova is much lower than the Galactic average, due to enhanced turbulence following previous 
novae for example, such a sustained source of cosmic rays will raise the local abundance. If 
efficient acceleration of particles to TeV energies in recurrent novae is commonplace, with 


spectral energy distribution harder than that of the Galactic cosmic-ray background x E ?^, 


the local contribution from novae would dominate at TeV energies over volumes > pc?. The 
size of the affected region will depend on the value of the diffusion coefficient, which can be 
constrained using measurements of the diffuse gamma-ray emission at energies ~ 10 — 100 
GeV (28). 

Our time-resolved gamma-ray emission measurements have implications for the origin of 
cosmic rays. Acceleration of cosmic rays to PeV energies in supernova remnants requires sub- 
stantial amplification of magnetic fields. Fast shocks (~ 10, 000 km s^!) propagating through 
the dense winds (M [Uwina ~ 10P kg m!) associated with the progenitors of supernova rem- 
nants from massive (2 8 Mo) stars, provide the only known environments where the required 
conditions can (in theory) be met (/7, 29). However, observational confirmation of this predic- 
tion has not been found. The detection of VHE gamma-rays from RS Oph provides an example 
of a Galactic accelerator reaching the theoretical limit for the maximum achievable particle 
energy via diffusive shock acceleration (77). If our results can be extrapolated to the most op- 
timistic supernova conditions, they support the prevailing model of Galactic PeV cosmic rays 


originating in supernova remnants from massive stars (77, 29). 
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Materials and methods 


H.E.S.S. Observations and Data Analysis 


HESS is an array of five Imaging Atmospheric Cherenkov Telescopes located in the Khomas 
Highland of Namibia, sensitive to VHE gamma rays in the energy range between a few tens of 
GeV and a few tens of TeV (30, 31). HESS consists of four 12 m-diameter Davies-Cotton 
telescopes (CT1-4), equipped with Cherenkov cameras (32) and have a field-of-view (FoV) of 
5°. The largest telescope (CTS) is located in the middle of the CT1-4 array, has a 28 m-diameter 
mirror and is equipped with a Cherenkov camera with 3.4° FoV (33). Observations of classical 
novae by H.E.S.S. are triggered by external alerts and follow-up observations are decided on a 
case-by-case basis taking into account if a transient signal is reported by Fermi-LAT, if the ejecta 
velocity is high (> 1500 km s7’), and if the peak optical magnitude is bright (my < 9). H.E.S.S. 
follows up on ~two nova candidates on average per year, subject to them satisfying at least one 
of the aforementioned conditions. Initial reports of an outburst of RS Oph on 21 — 9t? August 
2021 satisfied all three conditions, with a 6c detection by Fermi-LAT (34), ejecta velocity of 
> 2600kms~! (35, 36) and optical magnitude of my ~ 5.0, prompting H.E.S.S. observations 
with high priority. RS Oph was the first nova visible to H.E.S.S. to meet all three conditions at 
the same time. 

Follow-up observations of RS Oph with HESS. started on August 9*^, 18:17:40 (Coor- 
dinated Universal Time, UTC) and concluded on September 7*^, 19:47:31 (UTC). HESS, 
accumulated a total of 24.3 hours of observations over the first five nights after the explosion. 
After the full moon period, during which H.E.S.S. observations are paused, H.E.S.S. contin- 
ued RS Oph observations for a further 32.9 hours. The observations were conducted in the 
standard “wobble mode" where telescopes are pointed at alternate offsets of 0.5? from the po- 


sition of RS Oph (right ascension 17^50""13.165, declination -06?42/28.5" (J2000 equinox)). 


Night Tobs Livetime Significance Atmospheric transparency NSB noise level Telescopes 


(UTC) (hours) (c) 
09 Aug. 2021 18:17:40 3.2 5.8 (6.4) 0.90 1.0 CT1-5 
10 Aug. 2021 17:53:46 3.7 (2.8) 9.0 (7.1) 0.80 1.0 CT1-5 
11 Aug. 2021 17:44:08 3.7 9.8 (9.6) 0.65 1.0 CT1-5 
12 Aug. 2021 18:17:12 2.3 13.6 1.00 1.5 CT1-4 
13 Aug. 2021 17:44:43 2.8 10.5 (9.4) 1.10 2.5 CT1-5 
25 Aug. — 07 Sep. 2021 17:48:03; 19:47:31 14.6 (13.4) 3.3 (2.3) 0.96 1.0 CT1-5 


Table S1: Summary of H.E.S.S. observations and data set used in the analysis. The start time Tops and duration 
of observations is given for each night. Values are for the CT 1-4 stereo analysis, except numbers in brackets which 
denote CT5 mono analysis results. The atmospheric transparency is quantified using the Cherenkov transparency 
coefficient (37). The NSB noise level is given relative to the nominal Galactic level experienced by the telescopes. 
Only observations taken under good atmospheric conditions and low zenith angle < 35° were used in the data set 
from 25 Aug. 2021 to 07 Sep. 2021, leaving ~15 hours for analysis out of a total of 32.9 hours of observations. 
To achieve a low energy threshold, we only selected data for analysis with zenith angle < 35°, 
with a resulting average zenith angle of 22°. 

Throughout the first five nights, the data set is subject to strongly variable observing condi- 
tions. Due to increasing levels of moonlight, a spectral analysis of CTS data is performed for the 
first three nights only. On the second and third nights, atmospheric conditions were poor due to 
a higher-than-usual aerosol content in the atmosphere, whilst atmospheric conditions were good 
on the first, fourth and fifth night. We applied corrections to account for the night-to-night vari- 
ations in atmospheric conditions (see below). Furthermore, we conducted observations under 
low to moderate moonlight during nights 2 to 5, which increased the noise from night-sky- 
background (NSB) photons by a factor of 1.5 in night 4 and 2.5 in night 5 compared to dark sky 
conditions (c.f. Table[ST). 

We performed the analysis of the H.E.S.S. data using both events recorded by CT1-4 in 
stereoscopic (stereo) mode, and in a monoscopic (mono) analysis of CTS. The ring-background 
method (38) was employed for the extraction of skymaps for both analyses. The resulting sig- 
nificance maps are shown in Figure 1 for CT1-4. Background events were classified and re- 


jected using multivariate analysis techniques (39—41). The reconstruction of shower properties 


was performed using a pixel-based maximum likelihood technique (CT1-4 stereo, (42)) and a 


multivariate-based reconstruction (CT5 mono, (40, 41)). The spectral reconstruction and model 
fitting were performed with a likelihood-based procedure, using the GAMMAPY software pack- 
age (43-45) version 0.18.2. In the case of the CTS mono analysis, the reflected-region (38) 
method was used to define background control regions and to derive spectral results. 

Significant VHE gamma-ray emission was detected by H.E.S.S. from the direction of RS Oph 
on each of the five nights from 9'^ — 13*^ August 2021 (Table|S 1} in an online analysis and later 
confirmed offline by two independent analysis chains (39, 42, 46). The cross-check analysis 
employs an independent calibration, reconstruction and background suppression (46). On night 
4, CTS did not participate in observations for technical reasons, but it detected RS Oph in VHE 
gamma rays in the other four nights. 

RS Oph was observable again with H.E.S.S. under good low mooonlight conditions on Au- 
gust 25t” and observations continued as long as the source was visible to HESS until Septem- 
ber 2715. The analysis of these observations revealed a VHE gamma-ray signal from RS Oph at 
the 3 — 40-level in the main CT1-4 analysis as well as in the cross-check analysis. No sig- 
nificant emission was detected with the CT5 mono analysis in this late phase, consistent with 
its somewhat lower sensitivity. The corresponding flux upper limit is FuL( E > 250 GeV) = 
3.7x 107? cm~? s71. Table|S 1]summarizes the H.E.S.S. observations considered in our analysis 
below. 

The quality of the atmospheric conditions is quantified using the Cherenkov transparency 
coefficient (37) with lower values corresponding to lower transmission of Cherenkov light 
through the atmosphere, as reported in Table [ST] for the RS Oph observations. To correct for 
the different atmospheric conditions and resulting light-yield of the telescopes, the energy scale 
of the energy migration matrix as well as the effective areas have been scaled for individual 
observations with the atmospheric transparency coefficient for the CT1-4 stereo analysis. To 


validate this correction, the reconstructed Crab Nebula spectrum (after scaling the energy scale 
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Figure S1: Atmospheric transparency correction of Crab Nebula observations. Illustration of the effect of at- 
mospheric transparency correction to Crab Nebula observations taken during the same observing season as RS Oph 
for a range of atmospheric transparency conditions. The top panel illustrates the impact of the correction on the 
normalization œo of the reconstructed flux, ¢ = do - (E/1 TeV)-T. The bottom panel shows the reconstructed 
index I without and with correction. Vertical blue lines indicate the atmospheric transparency conditions corre- 
sponding to the five nights of RS Oph observations. 

of the instrument response functions (IRFs) with the transparency coefficient) is compared to the 
published H.E.S.S. spectrum, using the same IRFs as employed for the RS Oph observations. 
A good agreement between the spectral parameters of the corrected Crab Nebula spectrum and 
a previously published spectrum (30) is found for a range of atmospheric conditions, with the 
flux matching to within 10% (see Figure|S1}) 

For the CTS mono analysis, two different sets of IRFs tailored to the different atmospheric 
conditions have been used. A comparison study to the same Crab Nebula data set as described 
above was performed for the CT5 mono analysis, confirming that also here the Crab Neb- 
ula flux and gamma-ray spectral index is reconstructed correctly for the range of atmospheric 
transparencies experienced during the RS Oph observations (Figure|S 1). 

To test for spectral variation during the outburst, the H.E.S.S. data were analysed in nightly 


bins for the five nights with significant detection (see Table and fitted with a power-law 


model. Best-fitting parameters are given in Table [S2] for the first three nights using the CTS 


Data set oo Eo Index I 
[10-11 TeV-! cm-? s-1] [TeV] 
mono 09 Aug. 2021 14.9 + (2.7)stat. + (3-0) syst. 0.18 3.22 + (0.38)stat. + (0.20)syst. 
10 Aug. 2021 25.2 + (4.7)stat. + (5.0)syst. 0.18 4.01 + (0.48)stat. + (0.20) syst 
11 Aug. 2021 28.5 + (3.3)stat. + (5.7 )syst 0.18 — 3.15 + (0.23)sta + (0.20) syst. 
stereo 09 Aug. 2021 0.91 + (0.28) stat, + (0.14) syst. 0.35 4.24 + (0.75)stat. + (0.15) syst. 
10 Aug. 2021 1.90 + (0.32). + (0.38)sys 0.35 3.32 + (0.30) stat, + (0.15) syst. 
11 Aug. 2021 3.57 + (0.54)stat. + (0.54) syst. 0.35 4.08 + (0.4 2) stat (0.20) syst. 
12 Aug. 2021 3.00 + (0.33)stat + (0.45),u. 038 — 3.27 + (0.21) stat, + (0.15) syst 
13 Aug. 2021 1.77 + (0.25)sta + (0.35)sys 038 3.24 + (0.24) star + (0.15) syst 
stereo 25 Aug. 2021 - 07 Sep. 2021 | 0.238 + (0.080 )stat. + (0.036)sys. 0.35 3.33  (0.45)stat. + (0.15)syst. 


Table S2: Nightly spectral measurements of RS Oph. Best-fitting nightly spectral parameters from HESS, 


=P 
assuming a power-law model of the form du (&) , with derived systematic uncertainties. E is the reference 


energy, Do is the amplitude at the reference energy and T is the spectral index. 


mono analysis and for all five nights using the CT1-4 stereo analysis. Nightly spectral results 


are consistent between the two configurations. 
Systematic Uncertainties 


Multiple sources of systematic uncertainties contribute to affect the spectral measurements of 
the CT1-4 stereo and the CTS mono analyses. We summarise the various sources and their 
estimated influence on the reconstructed flux and gamma-ray spectral index. Systematic uncer- 
tainties stemming from Monte-Carlo extensive air shower hadronic interaction models, broken 
pixels of the Cherenkov cameras, and the live time of the data set are sub-dominant for both 
analyses, but have been taken into account (30). 

For the CT1-4 analysis, two different sets of cuts to select gamma-ray-like events have 
been applied to the night 1-5 data set, showing a systematic difference in the reconstructed flux 
normalisation and spectral index of ~5% and 0.1, respectively. The tests conducted using Crab 
Nebula observations under varying atmospheric conditions capture run-by-run differences in 
atmospheric conditions after correction, differences in the assumed and actual optical telescopes 
efficiency, and the scaling to correct for low atmospheric transparency itself. This results in a 
systematic error of the flux normalisation of 10%. The impact on the reconstructed spectral 


index varies from 0.05 (good, medium) to 0.15 (poor) atmospheric transparency. An additional 


10% systematic uncertainty of the flux normalisation is assumed for nights 2 and 5, when the 
Cherenkov cameras were operated at higher camera trigger settings and during moonlight. 

The uncertainty in the CTS mono analysis is derived from the following aspects. Following 
(47), an imperfect description of the background acceptance leads to a systematic uncertainty of 
15% on the flux normalisation and 0.15 on the spectral index for the CTS mono analysis. The 
Crab Nebula spectral comparison under comparable atmospheric conditions yields a systematic 
error on the flux normalisation of 5%, 10% and 15% for good, medium and poor atmospheric 
conditions, respectively. The uncertainty on the spectral index is at a level of 0.15 for the 
different transparencies. 

A comparison to the intrinsic spectrum of the Crab Nebula (44), shows that the overall 
energy scale uncertainty of the H.E.S.S. observations is at the 15% level. Table [S2] includes 
the resulting systematic errors of the RS Oph spectral measurements for nights 1 to 5 and the 
CT 1-4 stereo as well as the CTS mono analysis. The systematic error for the late follow-up 


observations is derived in the same way to that described above. 


Fermi-LAT Data Analysis 


We analyse Fermi-LAT data to investigate spectral evolution in bins of 24 hours, centred on 
20:00 UTC over the duration of the H.E.S.S. observations. We use Fermi-LAT (P8R3, pass 
8 release 3 (48)) data spanning from 9t* August 2021 to 7** September 2021, in the en- 
ergy range from 60 MeV to 500 GeV. We retrieve the data from a region of interest (ROI) 
defined by a radius of 15? around the position of RS Oph. The analysis of the LAT data de- 
scribed above is performed using the FERMIPY PYTHON package (version 0.18.0), based on the 
FERMI SCIENCE TOOLS (version 2.0.8) (49, 50). We analyzed SOURCE class events with a 
maximum zenith angle of 90? to eliminate Earth limb events. To derive the daily spectrum, the 


data are binned in 8 energy bins per decade and spatial bins of 0.1? size. (The re-binned points 
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in Figure 3 in the main article were determined only for plotting purposes, after the best-fitting 
parameters had been found.) The response of the LAT instrument is evaluated with the IRFs, 
version P8R3_SOURCE_V2) and we include in the model of the region all the LAT sources 
listed in the Fermi-LAT Fourth Source Catalog (4FGL, (5/)) in a radius of 20° around the po- 
sition of the nova. The contributions from Galactic and extra-galactic diffuse gamma-rays are 
described using the Galactic (gll_iem_v07) and isotropic (iso.P8R3_SOURCE_V2_v1) diffuse 
emission models. The models are available from the Fermi Science Support Center (FSSC) (52). 
The energy dispersion correction is applied to all individual models describing the diffuse and 
discrete gamma-ray emission in the ROI, except for the isotropic diffuse emission model. 

Spectral results obtained using 6 hour bins contemporaneous to the H.E.S.S. observations 
are found to be compatible with the results of the 24 hour bins, albeit with higher uncertainty. 
Additionally, we evaluate the systematic errors for the daily observations by applying the same 
analysis to a pre-flare 24 hour observation window. No significant gamma-ray signal is found at 
the position of the RS Oph nova in the pre-flare time interval, with a signal first detected during 
the optical rise time. This demonstrates that neither systematic effects in the Galactic diffuse 
emission model nor source confusion impact the Fermi-LAT flux measurement. 

To compute the light curve, the time range is split into a series of bins of varying duration. 
For the first five nights of the H.E.S.S. observations, contemporaneous 6h bins are used (and 
18 h bins in between these times). Outside of these five nights, we use bins of 24 h duration (or 
multiples thereof) as appropriate, given the statistics. To calculate the energy flux for a given 
time window, we model the sources” gamma-ray emission with a log-parabola spectral function 


—a-—flog(E/Eo) 
of the form ó(E) = do (#) SEH 


Eo , with parameters set to a spectral index a = 2.0, 


curvature? = 0.13 and reference energy Eu = 1 GeV, respectively. The flux normalisation ou 
is left free to vary to obtain the integrated flux in the 60 MeV to 500 GeV energy range. 


We use the derived Fermi-LAT spectral points for each observation period to perform joint 
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spectral fits with the H.E.S.S. mono and stereo data, using the GAMMAPY software package 
(version 0.18.2) (45). The different data sets were fitted simultaneously using a log-parabola 
spectral model; the best fitting parameters for each night are quoted in Table [53] The general 


trend is for the flux normalisation to decrease and the parabola to widen, over the course of the 


five nights. 
Data set do Eo [s B 
[1074 TeV-! cm~? s-1] [GeV] 
09 Aug. 2021 7.06 + 0.58 1.0 1.98 +0.04 0.19 + 0.01 
10 Aug. 2021 4.27 + 0.43 1.0 2.122 0.05 0.13 + 0.01 
11 Aug. 2021 3.69 + 0.42 1.0 2.01+0.06 0.13 + 0.01 
12 Aug. 2021 1.79 + 0.32 1.0 2.02+0.09 0.11 + 0.02 
13 Aug. 2021 1.94 + 0.32 1.0 2.05+0.07 0.12 + 0.02 


Table S3: Nightly log-parabola fits for RS Oph. Best-fitting spectral parameters from a joint fitting of the Fermi- 
LAT and HESS data with a log-parabola model. Eg is the reference energy, do is the amplitude at the reference 
energy, a is the spectral index and £ is the curvature. 


Supplementary Text 
Modelling of gamma-ray emission 


Recurrent novae are complex hydrodynamic processes accompanied by the formation of mul- 
tiple shock waves. The shocks are either internal, 1.e., formed by collision of different parts of 
ejecta, or external, i.e., formed in the circumbinary medium. Both internal and external shocks 
might be responsible for the production of gamma-ray emission. While internal shocks remain 
a hypothesis (see, however, (53)), external shock formation is unavoidable and the emission 
produced at the external shock is expected to be long-lasting. This was evident from the multi- 
wavelength observations of the 2006 RS Oph nova event (12-15, 54). We therefore focus on the 
external shock component. Similarly to supernova explosions, the external shocks of recurrent 
novae have different phases of evolution (55). Initially, when the circumbinary medium can- 
not substantially influence the explosion dynamics, the blast wave expands with approximately 


constant velocity, the so-called ejecta-dominated phase. Once a sufficient amount of material 
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is swept up, the shock decelerates, transitioning to the Sedov-Taylor phase. The shock decel- 
eration rate further increases in the radiative phase. For a point-like explosion in a radial wind 
with a mass density (p) profile p(r) x r~?, the shock radius (Rn) satisfies Ra œ 1?/? during 
the Sedov-Taylor phase (55). When cooling dominates the shock enters the radiative phase, 
with Ra, œ t'/?. Observations of the previous 2006 outburst of RS Oph indicate that after 
day ~ 6 the spectroscopic data is consistent with a shock speed usn x t 99 (12, 54). While 
usn X t~° over the first 3 weeks, it is inconsistent with the expected X-ray flux (13), indicating 
a self-similar solution does not apply. Follow-up imaging of the expanding shock from the 2006 
outburst (74) showed that the polar expansion continued at about 5 500 km s^! for at least 155 
days, apparently with minimal deceleration. This is consistent with 3D hydrodynamical simu- 
lations showing late-time expansion at > 4000 kms”! in the polar directions (7). This points to 
the potential complexity of the source geometry and difficulty in interpreting the measurements 
in the context of a model with a single shock velocity. The mid-latitude regions with slower 
shocks, denser gas and stronger radiative cooling, could dominate the spectral measurements. 
Assuming a similar evolutionary profile for the 2021 outburst, the source was detected with 
H.E.S.S. over the first five days primarily during the ejecta-dominated phase, during which the 
velocity is in the range (2 — 4) x 10? kms !. The velocity may have decreased by the time 
H.E.S.S. observations recommenced after the bright moonlight break. Optical spectroscopy ob- 
servations suggest that on day 2 the shock was faster than found in earlier measurements made 
on day 1 (71,36). This we attribute to the shock propagation in an inhomogeneous stellar wind, 


as found in previous 3D simulations (7). 


Shock conditions: We assume that the explosion origin is located at the white dwarf star, 
which moves around a red giant star on an orbit of radius forb © 1.48 au (7). The accretion onto 


the white dwarf proceeds through a dense disk in the orbital plane of the system. This results in 
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Figure S2: Schematic of external shock model. The explosion originates near the surface of the white dwarf 
(WD). Within one day, the shock is expanding as a bipolar blast wave moving orthogonal to the accretion disk, 
into the wind of the red giant (RG). The colour gradient of the shock indicates the expansion velocity, whilst the 
red gradient of the surrounding medium indicates the density of the RG wind. Material internal to the shock is 
shown in blue. Figure is not to scale. 

the quasi-spherical blast wave initially moving faster in directions perpendicular to the orbital 
plane, as two sections of a spherical shock, as illustrated in Figure[S2] This picture is consistent 
with observations of the previous nova event of RS Oph (/5). The period relevant for the 
H.E.S.S. observations likely covers the transition of the outflow expansion from the bipolar to 
quasi-spherical regimes (and also from the ejecta-dominated to the Sedov-Taylor phases), thus 


self-similar solutions (see, e.g. (56, 57)) cannot describe the entirety of this period consistently. 


We employ a simple numerical model that describes the motion of the ejecta (namely dis- 
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tance covered by the ejecta Rej, and its speed vej) and the dynamics of the explosion’s forward 
shock (its radius RA, and speed usn). These quantities are obtained under the following as- 
sumptions: (1) the density and pressure across the shocked region are constant, and determined 
by the Rankine-Hugoniot conditions (58); (ii) the upstream medium density is determined by 
the distance covered by the shock as 

z M 

Aro Dia + RA) 


Pup (S3) 


Here M is the mass-loss rate of the red giant star and vw is its wind speed, which depends on 


the distance as 


Bw 
R, 
Uw = Vo | 1 5 5 | (S4) 
V Vorb + Rin 


with və the terminal wind speed. We take the value R, ~ 0.35 au for the radius of the red 
giant (59) and adopt for the exponent of the wind acceleration D, = 3 (60). For the red giant star 
in RS Oph the wind parameter M /v,, is estimated to be in the range (5.7 — 7.4) x 10!!kgm ! 
(59). We therefore adopt a fiducial value of Min, = 6 x 10!! kg m™!. For a typical terminal 
wind velocity of v4 = (10 — 30) km s^! for red giants (61), this corresponds to a mass-loss rate 
of M z (1 — 3) x 1077 Ma yr-! (where M, is the solar mass). 

For a strong shock in single-atom polytropic gas, one can estimate the shocked gas (pas) 


mass-density as 


lau 


—2 
122 R2 
base 10” (re) gem ?, (S5) 


where we set vw = 20 km s71. The numerical value outside the parentheses in the shocked gas 
density expression equates to a proton number density of ng, & 5 x 10?cm ?. The correspond- 


ing pressure pas in the shocked region 


3 
Das © Tg Pastish (S6) 
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decelerates the ejecta according to 


duej 
Ti d 


= —4n Rina . (S7) 


Here mej is the ejecta mass, vej is ejecta speed (initially the ejecta is assumed to move with 


Vejo), and Raj is the distance traveled by the ejecta: 


t 


Sat) = f deal. (S8) 


0 


The radius of the forward shock HA is determined by the condition that the shocked region 


density is equal to the value given by Eq. (S5): 


Ve ps m (S9) 
where the total mass of the shocked region is given by 
Rat) 
Mas(t) = 4r f dr'r’? ouni zl). (S10) 


0 


We solve the above system of equations using an iterative method and obtain the time depen- 
dence of the forward shock speed, shocked region density, and ejecta speed. The kinematics of 
the forward shock and ejecta are shown in Fig. including the shock/ejecta speed (Fig. [S3A) 
and the path covered by the shock (Fig. [S3B). For the considered geometry the forward shock 
propagates in a medium with a constant density during the first ~ 12 hours, Fig. [S3B. After- 
wards, the upstream medium is characterized by a decreasing density, which results in a period 
of acceleration of the forward shock. The forward shock speed achieves its maximum value 
after approximately 3 days. The shock subsequently starts to decelerate, as expected from self- 
similar solutions in a medium with the density decreasing as pup «x 7? (56). 

For the shock speeds in Figure the post-shock temperature is > 10°K and the only 


effective radiative cooling process of the thermal plasma is bremsstrahlung. Fig. [S4]shows that 
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Figure S3: Temporal shock evolution. A: Model shock speed as a function of time since the explosion. Markers 
correspond to the expansion velocity estimated from the absorption center and the blue absorption edge of the Ha 
and HB P Cyg profiles obtained with Varese (36, 62) and SALT (//) spectroscopy of RS Oph. The shock speed 
is obtained in the framework of a single-zone model for initial ejecta velocity vejo = 3000 km s^. B: distance 
covered by the forward shock and distance between the center of the RG star and forward shock. 
the cooling time of the post-shock gas behind the forward shock (assuming the red giant wind is 
compressed by a factor of 4) is much longer than the expansion time of the nova for the duration 
of the H.E.S.S. observations, and so we expect the forward shock to be adiabatic, in contrast to 
models of slower internal shocks in classical novae (3). 

In the simplest treatment of diffusive shock acceleration at the external shock, assuming 


Bohm scaling (63), the acceleration rate for relativistic particles of charge q is determined by 


the shock speed and the post-shock mean magnetic field strength, Bsn as 


gp xem si» 
acc 

where 7) > 27 is a dimensionless parameter corresponding to the ratio of the relativistic particle- 

gyrofrequency to the scattering rate. Its inverse can be taken to characterize the acceleration 

efficiency and is a measure of the energy in the turbulent magnetic field. If non-linear magnetic 


field amplification occurs, 7 may fall well below 27. This should however be done with care, 


since magnetic field amplification is unlikely to occur on all scales simultaneously, which is 
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Figure S4: Cooling times. Bremsstrahlung cooling time (1.551) of the thermal plasma behind the forward shock as 
a function of the time since nova eruption (te) in our model (blue solid line), with the time elapsed since eruption 
also plotted for reference (orange dashed line). The pre-shock gas density, ng, is also shown (blue dotted line, 
right-hand y-axis). 

the essential principle of Bohm scaling (64). Particle acceleration to higher energies can be 
saturated either by escaping the acceleration region (confinement limited) or, in the case of 
electrons, when inverse Compton (IC) or synchrotron energy losses overtake the acceleration 
(cooling limited). The mean magnetic field strength at the forward shock is determined by the 
magnetic field of the red giant wind. We follow the Parker (65) description of a stellar wind 
expanding radially outwards from a rotating star, for which the magnetic field near the poles 
approaches a split-monopole (B œ r~?) and near the equator is swept into a Parker spiral with 
a toroidal-dominated field (B x r^!). Since the nova expansion is understood to be bipolar, 
and the rotation axis of the red giant should be similar to the orbital axis of the binary system, 


we adopt the split-monopole approximation for the magnetic field strength at distance R from 


the red giant: 
RA? 
BUR) = By (=) : (S12) 


where D, is the red giant surface magnetic field, typically of order 1 G (66). We consider two 
cases, adopting B, = 1 and 10 G (baseline and extreme values, respectively). The computed 


strength of the magnetic field in the blast-wave downstream is shown in Figure ER The wind 
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should be kinetic-energy dominated at a radius of a few times À,, so a limit on the magnetic 


field (assuming wind expansion at v,,) is 


M —1 
Be = 9G (=) (S13) 
R lau 


This is consistent with our assumed surface field of order (1 G) at R = R,. Substantially larger 


surface fields can be ruled out, as they would imply magnetically dominated winds at large 


distances. 
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Figure S5: Magnetic and photon field evolution. A: Model magnetic fields as a function of time since the 


explosion. B: Energy density in the optical photon field at the forward shock and in the magnetic fields as a 
function of time since the explosion. 


Maximum Energy: The VHE H.E.S.S. detection confirms the acceleration of particles to at 
least several TeV during the ejecta-dominated phase. The Hillas limit (67) for a particle of 


charge q — Ze in the above radial field (i.e. ignoring strong magnetic field amplification) is 


Ush Ush B, Fin = 
Ba eg EE 102 ( ) TeV S14 
d Ren B( Ben) 5000 km s-1 (3) (2) i ry 


assuming A = 0.35 au. In the first week the shock advances approximately 2 — 3 au per 


day, such that strong magnetic field amplification above that of the background value Bw at 
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the forward shock is essential to account for ongoing TeV particle acceleration > 3 days after 
the explosion. To facilitate acceleration, the magnetic fields are necessarily amplified upstream 
of the shock, which can only be triggered by energetic particle currents. For Bw x R~? the 
magnetic energy density falls rapidly, suggesting that only a small fraction of cosmic-ray energy 
needs to be converted to magnetic energy, for example via the Bell instability (68), to amplify 
the fields above background levels. This is an exact analogue of the core-collapse supernova 
scenario (77), for which there is a formalism to predict the maximum energy for a fast shock in 
a dense wind. We summarize that calculation (17) below. 

Neglecting shock curvature, and assuming isotropy of the particles, the differential acceler- 
ating flux (i.e. the rate at which particles increase their energy from E to E + dE) at a shock is 


determined by the shock compression (20) 


U(E) = SE Ne (S15) 


where Ng(E) is the differential non-thermal particle density, and Au = Ush — Uas, i.e. the 
difference between the upstream and downstream flow velocities. In the model (77), particles 
must generate their own self-confining waves at the expense of an escaping flux. Setting the 
escaping current density to the accelerating flux of the maximum energy particles when the 
shock was at a radius r i.e. ja«(r) = QU(Emax), and for simplicity assuming an Ng x E^? 


spectrum, it follows that 


Task 1 
m" ES Eese Sai (S16) 
where 
Au Us 1 
esc — S17 
GE us) 1og(Binax/ mc?) POR 


with Uer the energy density of cosmic-rays above the rest mass energy. The efficiency pa- 
rameter esc corresponds to the fraction of energy density flux processed by the shock that 


is lost to the upstream escaping energetic particles, and since Emax œ~ 107TeV, we expect 
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log(Emax/ Mc?) ~ 10. For acceleration efficiencies of U.,/ (E Pupu2,) zz 10%, approximately 
1% of the energy density flux processed by the shock is lost to the upstream. To determine if 
particles are confined, the fluctuations driven by the non-resonant instability (68) need sufficient 
time to grow to a level that can effectively scatter the highest energy particles. This requires 
many growth times of the fastest growing instability. The current that passes through a given 
fluid element at upstream radius À is the integral of the escaping current from the shock over 
the history of its expansion, diluted by a factor (r/ R)? (here we assume radial magnetic field 
and ignore deceleration of the shock). Using the definition of the maximum growth rate for 


cosmic-ray driven instabilities (68), and demanding at least 5 growth times, we find 


T maxdt = I J| — jedi > 5 (S18) 
PupC 


Using equation to determine Jesc, by the time the shock reaches a fluid cell located at radius 


R, the confinement condition is 


"TN ES = = a r-dr > 5R? (S19) 
Ro up Dosis 


where À, is the birth location of the shock. Differentiating both sides with respect to À, we find 


EQ, LES Em (820) 
20 Us C 


which is independent of both the magnetic field strength and shock radius. Eq. (S20) is the 


so-called confinement limit. This can also be used to estimate the effective scattering magnetic 
field Beg felt by the maximum energy particles found by equating the above expression for 


the confinement limit with the first equation in equation This results in a magnetisation 


B esc S 
GE (: | (S21) 
AT PupUs, 20 


parameter of 


This is only a representative value of the upstream field. Stronger magnetic fields concentrated 


on scales much less than the gyroradius of the maximum energy particles as measured in Beg 
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will not affect their acceleration. Further amplification can also occur in the post-shock medium, 
though this has only a minor effect on the acceleration rate, as the acceleration time is domi- 
nated by the upstream residence time of the diffusing particles. Because Emax depends only 
on the properties of the wind, the acceleration time for particles with energy Emax is always 
approximately equal to the instantaneous age of the shock. 

Using the shock and wind parameters discussed above for RS Oph, the maximum energy is 
several TeV. The prediction is consistent with the detection of TeV gamma-rays from RS Oph. 
The confinement limit is the dominant constraint for protons. Electrons must compete also with 
continuous IC, bremsstrahlung and synchrotron energy losses, so their maximum energy may 
be limited further. In the numerical leptonic model, the maximum electron energy does not 


exceed the confinement limit, but this still requires an escaping flux with €.,. œ 0.01. 


Radiative Losses: High-energy electrons can interact via different processes losing their en- 
ergy on time-scales comparable to the acceleration time. The typical cooling times for the 


synchrotron tsyn, IC (in the Thomson regime) fc 2. and bremsstrahlung tpr processes are 


E 1 / p\-? 
ts n x 4 aCe PEU ; 
SESCH (: xv) (5) g ee) 
E =l Woh s 
t ~ 15 | ——— ee : 2 
ee (; SH (; erg SH SS $929) 
and 
Nas =i 
fy, & 10° (=) . 24 
n o 10? cm-3 à SCH 


Here wp» is the energy density of target photons and na, is the number density of particles in the 
background plasma. These estimates show that for the conditions expected in the production 
region (see Fig. |S5), the IC mechanism plays the dominant role. For VHE electrons, Klein- 
Nishina suppression (69) decreases the rate of IC scattering, and synchrotron cooling might 


become dominant in this range. We therefore consider these two processes below. 
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As the radiative cooling becomes more important with increasing electron energy, we are 
interested in accurately estimating its contribution for VHE electrons. TeV electrons efficiently 
up-scatter target photons with energy below the limit set by the Klein-Nishina cutoff, < 0.25 eV. 
Thus, optical and infra-red (IR) observations provide information about the intensity of the soft 
photon fields, which serve as a target for IC losses and emission. Assuming these photon 
fields result from thermal free-free emission, the expected spectrum in this energy band can be 
approximated as 


Sy ocu * BLUT (0), (S25) 


where B, is the Planck function, T is the (time-dependent) temperature of the gas producing 
free-free emission, and the factor v~* appears due to free-free opacity. For the typical post- 
shock temperatures T > 108 K, we are interested in the lower energy part of the spectrum, 
hv < kpT, and the blackbody spectrum can be approximated by the Rayleigh-Jeans law, 
B, x v?. Here h and kg are Planck and Boltzmann constants, respectively. The relevant part 
of the spectrum should have a flat spectral energy distribution (SED). In our model calculations 
we directly use the time-dependent information from the optical B, V, R, and I bands (70), and 
below the I band we approximate the target photon spectrum by a flat SED component. This 
approximation agrees with the IR spectra measured from RS Oph after its eruption in 2006 (70). 
The time-dependence of the photon target energy density is shown in Figure[S5B. 


The evolution of the non-thermal particles follows the acceleration-losses equation: 


dE dE 
dt dt 


| dE 
dt 


; (S26) 


losses 


acc 


where the acceleration rate is given by equation (S11) and for electrons the loss term ac- 


— dE dE . . . 
= T alie Relativistic 


counts for both synchrotron and IC mechanisms, als 


d 
dt |losses 
bremsstrahlung is important at early times, but negligible compared to synchrotron and IC by 


the time the H.E.S.S. observations commenced. Similarly, for the expected densities, the energy 
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Figure S6: Particle energy and luminosity evolution. A: Maximum energy of accelerated electrons and protons 
at the forward shock. B: Energy crossing the forward shock per time unit, La, and observed optical-IR luminosity 
of the nova, Lon, The dip occurring around day 1 was in all optical bands. 
losses for protons are negligibly small. The maximum energy determined from the solution to 
equation for the two cases considered is shown in Figure [S6] For this plot we adopted 
n = 2n both for electrons (ne) and protons (np). In the case of protons, we assume the maximum 
energy follows the confinement limit described above with £4, = 107?. For protons the max- 
imum energy does not depend on the magnetic field strength (except for a short time interval 
after the onset of the acceleration, which is not distinguishable on the scale of the figure). The 
maximum energy grows rapidly at early time due to the strong magnetic fields close to the star, 
and the 7 = 27 example shown in Figure [S6]likely over-estimates the acceleration rate. Follow- 
ing the rapid growth, the maximum energy saturates at the confinement limit and follows this 
value for later times. Thus non-linear magnetic field amplification is implicit in the numerical 
model. The transition from fast to slow cooling is also evident in the electron spectrum . 

The obtained maximum energy determines the spectral cut-off of the non-thermal particles, 


which is injected downstream at each moment: 


Be/p 
dau ` J AE /rexp |- (z " ) | | E > Enine/p 


= max,e/p 
B 0 Den 


(S27) 
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Here Emax,e/p is obtained by solving Eq. subject to the limit in Eq. (S20), for electrons 
and protons; Emin,e/p is a free parameter set to 100 MeV for leptons and 2 GeV for protons to 
avoid non-relativistic effects. The normalization coefficient À is computed by the condition 
that a fixed fraction, k./,, of the energy flux through the forward shock, Lsn, is transferred to 


non-thermal particles: 


p= Ke/pLsh 
e/p ` oo S 
f | dEE'-* exp |- | 


min,e/p 


(S28) 


E 


max,e/p 


E 
The energy flux through the shock is determined by the upstream medium density and the 
forward shock speed: 


La = 2T R? pap, (S29) 


The model dependence of this parameter is shown in Figure S6B. 
The spectrum of non-thermal particles, n = dN/dE, is described by the single-zone evolu- 


tion equation 


On o : dNinj 
Ot + dE |n — dEdt , (S30) 


where the injection term is given by Eq. and the initial condition is n|, , = 0. 

The non-thermal particle spectra expected in the source at days 1, 2, 3, 4, and 5 after the 
explosion are shown in Figure In the case of non-thermal protons the time evolution of 
the spectrum is simple. The impact of cooling on the proton spectrum is minor and we expect 
an accumulation of non-thermal protons accompanied by minor spectral evolution. In contrast, 
the development of the non-thermal electrons is quite complex, and the spectral properties are 
affected by the strength of the magnetic field. There are features caused by the decrease in 
energy loss rates with radius: increasing maximum energy and spectral hardening caused by 
the transition from the fast cooling to slow cooling regimes. For the chosen parameters, we do 


not expect acceleration of electrons to very-high energies during the first day after the explosion 
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Figure S7: Single-zone particle spectra evolution. A: Electron spectra for five instants in time post-outburst and 
two different values for the strength of the red giant surface magnetic field. B: Spectra of non-thermal protons for 
five instants in time (the impact of magnetic field strength is minor, due to confinement limited saturation). 

for the case of weak external magnetic field. On the other hand, the flux level in the VHE band 
is smaller than that at GeV energies thus a sufficient VHE flux can be produced by particles with 
energy exceeding the cutoff energy, in the tail of the distribution. Therefore for the radiation 
calculations we consider the case with B, — 1 G, which as discussed above is more consistent 
with the value expected at the RG surface. 

In Figure|S8] we compare the model IC spectra produced by the shock-accelerated electrons 
to observational data obtained with Fermi-LAT in the GeV band and H.E.S.S. data obtained in 
the VHE regime. We set &, — 0.1 and compute the IC spectra for five different instants in time 
post-outburst: t = 1, 2, 3, 4, and 5 days after the explosion. The IC model under-predicts the 
flux on night 1; although an increased IC flux for the first night can be obtained by increasing the 
density of the wind in the model, this change would also accelerate the evolution, with the wind 
deceleration commencing between days 1 and 2 post explosion. A denser wind is therefore 
problematic due to the rapid deceleration of the shock (incompatible with multi-wavelength 
observations) and the need to introduce a suppression of the p-p channel to maintain dominant 


IC (inconsistent with the theory of shock acceleration). 
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In Figure[S9]we show the result of hadronic emission spectra (77) due to collisions between 
the non-thermal protons and shock compressed external medium for the same times, using 
Kp = 0.5. Here we assume that the gamma-ray emission is produced in a region close to the 
external shock where the density is twice higher than the averaged value (1.e., the filling factor 
is 50% ). The inverse of the acceleration-rate efficiency for electrons was set to 7. = 107 and 


the maximum energy of protons was calculated for £4. = 107? and n, = 307. 
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Figure S8: Nightly Inverse Compton spectra. A-E IC emission of shock-accelerated electrons for nights 1-5 
respectively, together with spectral measurements obtained with H.E.S.S. (filled points) and Fermi-LAT (hollow 
points). The model parameters are summarized in Table[S4] 

The local intense photon fields, especially during the first days after the explosion, require 
consideration of y — y attenuation. We provide some general estimates, which allow us to define 


the photon frequencies and epochs when ^ — y absorption might be important. For a source 


with luminosity L and size Rx we can put the following upper limit for the optical depth, 7, for 
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Figure S9: Nightly hadronic emission spectra. Same as Fig. but for proton-proton emission. The model 
parameters are summarized in Table[S4] 


photons with energy e: 
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Here 0.207 is the maximum cross-section for the pair creation process expressed in the Thom- 


son cross-section value (or). This maximum is achieved for chijarget ~ 3.5m2c*, where hiarget 
is the target photon energy. The attenuation of TeV gamma-rays occurs due to the interaction 
with optical — IR photons, and X-rays provide the dominant target for GeV photons. 

Equation (S31) shows that the y — y attenuation is the most relevant for VHE photons. 
For this gamma-ray energy range, optical and IR photons provides the dominant target, and 
below we compute the corresponding optical target photons. During the first several days after 


the eruption the typical optical — IR luminosity remains at the level of (1 — 5) x 10?"ergs ^! 
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(see Fig. [S6]and (10)), while the typical X-ray luminosity is fainter. For example, X-ray fluxes 
at the level of 1 keV cm ?s ^! were reported at days 3 — 6 after the eruption in 2006 (13). 
For the adopted source distance of 1.4 kpc this corresponds to the X-ray luminosity of ~ 3 x 
109 ergs-!. This implies that the y — ~y attenuation might be important for GeV photons only 
during the first hours after the eruption, when the source size is small, Ra, << 1 au. 

To estimate its impact we compute also the attenuation factor, exp [—7], averaged over the 
forward shock. The resulting factors, again assuming a distance of 1.4 kpc are shown in Fig- 
ure[S10]for five different epochs: t = 1, 2, 3, 4, and 5 days after the explosion. The absorption 
is not sufficient to account for the observed difference between the HE and VHE bands, and 
therefore must reflect a feature in the particle population responsible for the gamma rays. The 
impact of the y — y attenuation is shown in Fig. where it can be seen that the spectrum 
transformation by the absorption is minor. We note that the spectral energy distributions in 
Figs. and [S9] account for the attenuation factor obtained for the source distance of 1.4 kpc 
shown in Fig. [STOA. 

The idealised single zone model presented here does not accurately account for the inho- 
mogeneous structure of the post-shock medium, in which strong magnetic field and gas density 
gradients naturally develop. The purpose of this model is therefore not to precisely fit the data, 
but rather to reproduce essential features of the system’s evolution, characterised by a minimally 


sufficient number of free variables. 
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Figure S10: Effect of gamma-gamma attenuation. A: Attenuation factor due to the gamma-gamma absorption 
on the thermal photons produced by the explosion assuming a distance of 1.4 kpc. B: The impact of the gamma- 
gamma absorption on the emission spectrum from night 1 for the source assumed distance of 1.4kpc. The data 
points are the same as in Fig.[S9] 


Parameter Symbol, unit p-p model IC model 
Acceleration slope electrons Qe - 2.2 
Acceleration slope protons Qp 2.2 - 
Cutoff exponent electrons Be - 0.5 
Cutoff exponent protons Bp 0.5 - 
Fraction of energy in electrons Ke 0 3% 
Fraction of energy in protons Kp 50% 0 
Acceleration efficiency of electrons Ne - 107 
Acceleration efficiency of protons Np 307 - 
Escape efficiency bese 10-2 - 
Electron low energy cutoff Emin. - 102m.c? 
Proton low energy cutoff Pain 2.mpc? - 
RG surface magnetic field Bs, G I. 

RG radius, au Rx, au 0.35 

RG mass-loss rate M/ow, g cm! 6.3 x 1011 

WD orbit radius Torp, au 1.48 
Distance from Earth kpc 1.4 

Ejecta initial speed Dein, kms” T 3000 

Ejecta mass Mej 10-7 Mo 


Table S4: Numerical model parameters. Assumed parameters of the model used in equations through (S30) 
above, for the two cases of proton-proton (p-p) emission (see Figure [S9) and a leptonic Inverse Compton (IC) 
model (see Figure[S8) respectively. 
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